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ABSTRACT 

We study the early dynamical evolution of young, dense star clusters using Monte Carlo simulations 
for systems with up to iV = 10^ stars. Rapid mass segregation of massive main-sequence stars and the 
development of the Spitzer instability can drive these systems to core collapse in a small fraction of 
the initial half-mass relaxation time. If the core collapse time is less than the lifetime of the massive 
stars, all stars in the collapsing core may then undergo a runaway collision process leading to the 
formation of a massive black hole. Here we study in detail the first step in this process, up to the 
occurrence of core collapse. We have performed about 100 simulations for clusters with a wide variety 
of initial conditions, varying systematically the cluster density profile, stellar IMF, and number of 
stars. We also considered the effects of initial mass segregation and stellar evolution mass loss. Our 
results show that, for clusters with a moderate initial central concentration and any realistic IMF, the 
ratio of core collapse time to initial half-mass relaxation time is typically ~ 0.1, in agreement with 
the value previously found by direct iV-body simulations for much smaller systems. Models with even 
higher central concentration initially, or with initial mass segregation (from star formation) have even 
shorter core-collapse times. Remarkably, we find that, for all realistic initial conditions, the mass of 
the collapsing core is always close to ~ 10~^ of the total cluster mass, very similar to the observed 
correlation between central black hole mass and total cluster mass in a variety of environments. We 
discuss the implications of our results for the formation of intermediate-mass black holes in globular 
clusters and super star clusters, ultraluminous X-ray sources, and seed black holes in proto-galactic 
nuclei. 

Subject headings: Black Hole Physics — Galaxies: Nuclei — Galaxies: Starburst — Galaxies: Star 
Clusters — Methods: N-Body Simulations, Stellar Dynamics 



1. INTRODUCTION 
1.1. Astrophysical Motivation 

It is now well established that the centers of most 
galaxies host supermassive black h oles (BH) with masses 
in the range ^ 10^ - 10^ Mq fF errarese et"al]l2001b 
Eormcndv & Gcbhardt 2001). The evidence is partic- 
ularly compelling for a BH of m ass c± 4 x 10^ Mfr, 
at the center of our own Galaxv llEckart et all 120021: 
iGhez et al.ll200(l 1200^ iSchodel et al.ll2002() . Dvnamical 
estimates indicate that, across a wide range, the cen- 
tral BH mass is about .1% of th e spheroidal compo- 
nent of the host galaxy l)Holll998|) . A related correla- 
tion may exist with the total gravitational mass of the 
host galaxy (ba sically the mass of its dark matter halo) 
iFerrares^ (|2002l) . An even tighter correlation is observed 
between the central velocity dispersio n and the cen- 
tral BH mass (Fcrrarcsc fc Merritljl20'ont iGebhardt et all 
I2OOO:; .Tremaine et al 2002.) . 

Theoretical arguments and recent observations suggest 
that a cent ral BH may also exist in some globular clus- 
ters (jvan d cr Marcl 2001, 2003). In particular, recent 
HST observat i ons and dy namical modeling of Ml 5 by 
IGerssen et all (|2002l 1200.'^ yielded results that are con- 
sistent with the p resence of a central ma ssive BH in this 
cluster. Similarlv lGebhardt et al.l l)2002() have argued for 
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the existence of an even more massive BH at the center 
of the globular cl uster Gl in M31. H owever. TV-body 
simulations (Baun igardt et al.ll2003albD suggest that the 
observations of M15 and Gl could be explained equally 
well by the presence of many comp act objec ts near the 
center without a massive BH l|van deT Mare]|l200l|) . 

When the correlation between the mass of the cen- 
tral BH and the spheroidal component in galaxies is ex- 
trapolated to smaller stellar systems like globular clus- 
ters, the inferred BH masses are ~ 10'^ — 10^ Mq, much 
larger than a ^ 10 Mq stellar-mass BH, but much smaller 
than the ^ 10^ — 10^ Mq of supermassive BH. Hence, 
these are called intermediate-mass black holes (IMBH). 
If some globular clusters do host a central IMBH the 
question arises o f how these objects were formed (for re - 
cent reviews see lR,asio et al.ll2Cl03t Ivan der Ma rel"2003') . 
One natural path for their formation in any young stel- 
lar system with a high enough density is through run- 
away collisions and mergers of massive stars follow- 
ing core collapse. These runaways could easily oc- 
cur in a variety of observed young star clusters such 
as the "young populous clusters" like the Arches and 
Quintuplet clusters in our Galactic center and the "su- 
per star clusters" observed in all st arburst and galac- 
tic merger environments (see, e.g.. iFiger et al.l Il999al: 
'Galla gher fc Sm ith 1999). The Pistol Star in the Quin- 
tuplet cluster (Figcr ct al. 1998} may be the product 
of such a runaway, as d emonstrated recently by di- 
rect iV-body simulations l)Portegies Zwart fc McMillanI 
l2002t) . A similar process may be responsible for the for- 
mation of seed BH in proto-galactic nuclei, which could 
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then grow by gas accretion or by mer ging with other 
IMB H formed in young star clusters ijEbisuzaki et al.l 
I200H iHansen fc Milosavlievicll200^ . Further observa- 
tional evidence for IMBH in dense star clusters comes 
from recent Chandra and XMM-Newton observations 
of ultra-luminous X-ray sources, which are often (al- 
though not always) associated with young star clus- 
ters and whose high X-ray luminosities in many cases 
suggest a compa ct object mass of at leas t ^ 10^ M^t^ 
jEbisuzak i et all 120011: iKaaret et all I200H |Miller et aL| 
12003) , although beamed emission by an accre ting stellar- 
mass BH provides an alternati ve explanation l)King et al.l 
l200TllZeza.s Fahhia,nnll200l. 

When they are born, star clusters are expected to con- 
tain many young stars with a wide range of masses, 
from ~ O.IM© to ~ lOOM©, distributed according to a 
Salpeter-like initial mass function ( IMF) llClarkell200.^) . 
Inspired by an early version from iReed l)1984|) for the 
formation of supermassive BH, we show in Figure ^ a 
diagram illustrating the two main scenarios leading to 
the formation of an IMBH at the center of a dense star 
cluster. The early stages of dynamical evolution are 
dominated by the stars in the upper part of the mass 
spectrum. Through dynamical friction, these heavy stars 
tend to concentrate toward the center and drive the sys- 
tem to core collapse. Successive collisions and mergers of 
the massive stars during core collapse can then lead to a 
runaway process and the rapid formation of a very mas- 
sive object containing the entire mass of the collapsing 
cluster core. Although the fate of such a massive merger 
remnant is rather uncertain, direct "monolithic" collapse 
to a BH with no or little mass loss is a lik ely outcome, 
at le ast for sufficiently small metallicities ijHeger et al.l 
|2QQ2)- A_n essential condition for this runaway to oc- 
cur is that the core collapse must occur before the most 
massive stars born i n the cluster end their lives in su- 
pernova explosio ns ijPortegies Zwart & McMillai] 120021: 
iRasio et al.il200a) . 

The accumulation at the center of a galaxy of many 
IMBH produced through this runaway process in nearby 
young star clusters (like the Arches and Quintuplet clus- 
ters in our Galaxy) provides an interesting new way 
of building up the mass of a centr al supermassive BH 
l|Port.egies Zwart fc McMillanll2003) . It is possible that 
this process of accumulation is still ongoing in our own 
Galaxy (Hansen & Milosavlicvic 2003). These ideas have 
potentially important implications for the study of su- 
permassive BH by the Laser Interferometer Space An- 
tenna (LISA), since the inspiral of IMBH into a super- 
massive BH provides the best source of low-frequency 
grav itational waves for dire ct study of strong field grav- 
ity ifCutler fc Thorn e 2002). 

In the alternative scenario where massive stars evolve 
and produce supernovae before the cluster goes into core 
collapse, a subsystem of stellar-mass BH will be formed 
(Figure . As demonstrated in Section 14.51 the mass 
loss from supernovae provides significant indirect heat- 
ing of the cluster core, delaying the onset of core col- 
lapse until much later, after the stellar remnants un- 
dergo mass segregation. The final fate of a cluster with 
a component of stellar-mass BH remains highly uncer- 
tain. This is because realistic dynamical simulations 
for such clusters (containing a large number of black 
holes and ordinary stars with a realistic mass spec- 



trum) have yet to be performed. For old and relatively 
small systems (such as Galactic globular clusters), com- 
plete evaporation is likely (with all the stellar-mass BH 
ejected from the cluster through 3-body and 4-body in- 
teractions in the dense core). This is expected theo- 
retically on the basis of simple qualitative argum ents 
ijKulkarni et alJIlQQStlSigurdsson fc Hernauistill993|) and 
has been demonstrated recently by direct TV-body simu- 
lations for very small systems cont aining only ^ 10 BH 
ijPortcgics Zwart fc McM illan 20001. However, for larger 
systems (more massive globular clusters or proto-galactic 
nuclei), contraction of the cluster to a highly relativis- 
tic state could again lead to successive mergers (driven 
by gravitatio nal radiation) and the formation of a singl e 
massive BH (iLeeHlQQ.'il mM IQuinlan fc Shapird HQSt?). 
Moreover, it has been recently suggested that, if stellar- 
mass BH are formed with a broad mass spect rum (a likely 
outco me for stars of very low metallicity; see lHeger et alJ 
l2003j) . the most massive BH could resist ejection, even 
in a system with low escape velocity such as a globu- 
lar cluster. These more massive BH could then grow by 
repeatedly forming binaries (through exchange interac- 
tions )_^Tth_othe£j5H and merging with their compan- 
ions ijMiller fc Ham ilton 2002). However, as most inter- 
actions will probably result in the ejection of one of the 
lighter BH, it is unclear whether any object could grow 
substantially through this mechanism before running out 
of companions to merge with. 

1.2. Core Collapse and the Spitzer Instability 

The physics of gravothermal contraction and core col- 
lapse is by now very well understood for single- component 
systems (containing all equal- mass stars). In particular, 
the dynamical evolution of an isolated, single-component 
Plummer sphere to core collapse has been studied ex- 
tensively and has become a testbed for all numerical 
code s used to comp ute the evolution of dense star clus- 
ters l|Aarseth et al.lll974tl(^ersz fc Sour zemll 19941) . This 
evolution can be visualized easily using Lagrange radii, 
enclosing a fixed fraction of the total mass of the sys- 
tem (see, e.g., Fig. 3 of Joshi ct al. 2000 and Fig. 5 of 
iFreitag fc Benzll2001|) . As the system evolves, the in- 
ner Lagrange radii contract while the outer ones expand. 
This so-called gravothermal contraction results from the 
negative heat capacity that is a common property of all 
gravitationally bound systems (Elson et al. 1987). In the 
absence of an energy source in the cluster core, and in the 
point-mass limit (i.e., neglecting physical collisions be- 
tween stars) , the contraction of the cluster core becomes 
self-similar and continues indefinitely. This phenomenon 
is known as core collapse and its universality is very well 
established (see, e.g.. IFreitag fc Benzil2d0ll Table 1 and 
references therein). 

When the system contains stars with a variety of 
masses, the evolution to core collapse is accelerated. 
This has been demonstrat ed by the earliest A^-body and 
Monte Carlo simulations l)Aarsethlll966t lHenonlll971at 
lWielenlll97,'iD . To illustrate this behavior, and as a pre- 
view of the results presented later in Sec. 4, we show 
in Figure |21 the evolution of a cluster described initially 
by a simple Plummer model containing 1.25 x 10^ stars 
with a broad Salpeter IMF between mmin = 0.2 Mq 
and mmax = 120 M©. Core collapse occurs after < 
O.liih(O), where trh(O) is the initial half-mass relaxation 
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Fig. 1. — Various possible scenarios for the early dynamical evolution of a dense star cluster with a realistic IMF. The path for the 
formation of an IMBH through core collapse and runaway collisions (studied in this paper) is indicated by a thicker line. An alternative 
scenario involves successive mergers of stellar-mass BH binaries driven by a combination of dynamical interactions and gravitational 
radiation (left side of the diagram). For a runaway to occur, the core collapse time tec must be smaller than the stellar lifetime 4* of the 
most massive stars in the cluster. High-velocity disruptive collisions, the formation of a very extended and diffuse merger remnant, or the 
accumulation in the cluster of gas released by stellar winds and supernova explosions could lead to the formation of a complex system 
containing stars embedded in a dense gas clouds. The final fate of such a system is highly uncertain. 



time (see eq. © below). In sharp contrast, core col- 
lapse in a single-component Plummer model occurs after 
> 10ti.h(0) (a well-known result). Thus the presence of a 
broad IMF can dramatically accelerate the evolution of 
the cluster to core collapse. 

This acceleration of the evolution to core collapse is 
due to the changing nature of energy transfer in the 
presence of a wide mass spectrum. Relaxation pro- 
cesses tend to establish energy equipartition (see, e.g., 



iBinnev fc Tremain3ll987l Sec. 8.4). In a cluster where 
the masses of the stars are nearly equal, this can be (very 
nearly) achieved. The core collapse is then a result of en- 
ergy transfer from the inner to the outer pa rts of the clus- 
ter, leading to gravotherrn al contraction ijLarso 
iLvnden-Bell fc Woodlll96^ . A large difference between 
the masses of the stars allows a more efficient mechanism 
for energy transfer. In this case, energy equipartition 
would tend to bring the heavier stars to lower speeds. 
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Fig. 2. — Evolution of single-component and Salpeter-IMF 
Plummer models to core collapse. Lagrange radii enclosing a con- 
stant mass fraction, indicated on the left, are shown as a function 
of time. The radii are in units of the initial half-mass radius of 
the cluster and time is in units of the initial half-mass relaxation 
time. The solid lines are from our Monte Carlo simulation of a 
Plummer model containing 1.25 X 10® stars with a Salpeter IMF 
(within TTlniin — O.2 7Vf0 and TTlmax — 120 M0). The dotted lines 
show the result of a gaseous model simulation of the same cluster, 
with 50 discrete mass components approximating a Salpeter IMF 
(within the same mass limits). See Sec. 2 for more discussions of 
the various methods. For comparison, the dashed lines show the 
evolution of a single-component Plummer Model (computed with 
our Monte Carlo code). Note our key result: the core collapse time 
is more than two orders of magnitude shorter for a cluster with a 
realistic IMF. 



However, as a result, the heavier stars sink to the center, 
where they tend to gain kinetic energy, while the lighter 
stars move to the outer halo. This process is called "mass 
segregation." As the mass segregation proceeds, the core 
contracts and gets denser, leading to a shorter relaxation 
time, which in turn increases the rate of energy transfer 
from heavier to lighter stars. In typical cases this evolu- 
tion eventuall y makes the h eavier stars evolve away from 
equipartition il SDitzeil lT969>). 

The fundamental inability of the heavier stars to es- 
tablish energy equipartition with the lighter stars in a 
system with a continuous mass spectrum is similar to the 
Spitzer " mass-segregatio n instability" in two-component 
clusters. 'Soitzer' f 19691), using analytic methods and a 
number of simplifying assumptions, determined a simple 
criterion for a two-component system to achieve energy 
equipartition in equilibrium. If the mass of the lighter 
(heavier) stars is mi {1712) and the total mass in the hght 
(heavy) component is A/i (M2), then Spitzer's criterion 
can be written 



3/2 



<0.16. (1) 
iWatters et all (|2000j) . using numerical simulations, ob- 



tained a more accurate empirical condition, 

2.4 

A = ( ^ 1 ( — 1 < 0.32. 
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When this stability criterion is not satisfied, energy 
equipartition cannot be es tablished between heavy and 
light stars. .Spitzerl ljl969|) noted that the equililDrium 
would not be achieved for realistic mass spectra, because 
there is always sufficient mass in high-mass stars that, 
through mass segregation, they can form a subsystem 
(near the center of the cluster) that decouples dynami- 
cally from the lower-mass stars. This was later supported 
by rn ore detailed theoretical st udies ([I nagg^ fc Saslawl 
ITmiSaslaw He Younsll1971blVishnia,Hn97SD . 

Here we carry out a simple calculation to show that 
a typical system, with a Salpeter IMF (see Sec. 3.2) ex- 
tending from 0.2 Mq to 120 Mq, when viewed as a two- 
component cluster of "light stars" and "heavy stars," 
does not satisfy any of the simple stability criteria. Let us 
separate the stars into these two components according 
to some arbitrary boundary: we group all stars lighter 
than mb in the first component, and all stars heavier 
than mb in the second component. We use mi and Mi 
to denote the values of the average and total mass in the 
first component, and m2 and M2 similarly for the second 
component. We obtain, with mb in solar mass. 
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When these values are used in equations Q and ij^Jl we 
see (Figure that the stability criteria are almost never 
satisfied for any value of mb , except when it is extremely 
close to the maximum mass (as expected, since the model 
reduces artificially to a single-component system in this 
li mit, with all sta rs having the average mass). 

iVishniad l)1978l) . devised a criterion genuinely adapted 
to clusters with a continuous mass spectrum, under the 
ad-hoc assumption that the shape of the density distri- 
bution does not depend on the stellar mass. He derives 
the following necessary condition for stability: 
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M2 
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< 1, with /3 ~ 0.5, (5) 



for all mb (his eq. 17). Figure 13 shows that this condition 
cannot be satisfied either3. These results suggest that for 
any Salpeter-like IMF, one can always find a collection of 
stars that have large enough average and total mass that 
they will dynamically decouple from the lighter stars in 
the cluster. Once these decouple they form a subsystem 
with a much shorter relaxation time and consequently 
evolve to core collapse very rapidly. 

The half-mass relaxation time (relaxation time at the 
half-mass radius rh) for a cluster of N stars is given by 

^ The applicability of Vishniac's cri terion appears somewhat 
questionable. Through FP simulations. Ilnagaki &: Saslawl 119851) 
have shown that an IMF exponent a > 6.0 fsee~Sec~ l?2t is re- 
quired for central equipartition to set in before core collapse while 
Vishniac's criterion predicts that a > 3.5 is sufficient. 
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Fig. 3. — Various criteria for the realization of equipartition in 
a cluster with a Salpeter IMF. According to these, energy equipar- 
tition in dynamical equilibrium can only be achieved if the cor- 
responding curve lies entirely under the dotted horizontal line. 
Conditions by Spitzer and Watters et al. were devised for two- 
component models, so we divide the mass function into stars lighter 
and heavier than some arbitrary boundary value ra^ and we evalu- 
ate the stability conditions for any value of rrit,. Vishniac's analysis 
is genuinely adapted to a continuous mass function. He derives a 
necessary condition for stable equipartition that must be obeyed 
for any stellar mass mj, in the range covered by the mass spectrum. 
See text for details. 



l)Spitzedll987l eq. 2.63): 
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where p is the mass density and 7c 0.01 in the Coulomb 
logarithm. Let us assume that a coUapsing subsystem 
is formed by stars that constitute 1% of the total mass, 
and that aU of them come from the uppermost part of the 
mass spectrum. For a Salpeter IMF with mmin = 0.2 A/q 
and mmax = 120 Mq, the number of stars in this subsys- 
tem is then A^sub ^5 10^^ A^. At the time of dynamical 
decoupling the central density of the subsystem must be 
comparable to the central density of the overall cluster. 
Therefore we conclude that the relaxation time is around 
three orders of magnitude smaller for the collapsing sub- 
system of heavy stars at the onset of instability. So, 
essentially, the heavy stars will go into core collapse as 
soon as they start dominating the mass density near the 
center. 

If the cluster starts its evolution with heavy stars 
distributed throughout, the timescale for core col- 
lapse will be determined by their mass segregation 
timescale. The process of mass segregation for the heav- 
iest stars in the cluster is driven by dynamical friction 
l|Binnev fc Tremaindll987l Section 7.1). The timescale 
for mass segregation and the onset of core collapse there- 
fore depends primarily on the mass ratio between the 
dominant heavy stars and the lighter bac kground stars 
l|Chernoff Weinberg! ITflflil ISnitzeHrmfiflf^ . For simple 
dynamical friction in an effective two-component model, 
one would then expect the core collapse time to be com- 



parable to tDF = {'m))trh- Here mb (-C ^Mmax, 

the upper mass limit of the IMF) should be the mass 
above which the IMF contains a large enough number of 
massive stars to define a "collapsing subsystem." If this 
number is ~ lO'^ and the total number of stars N ~ 10^ 
we get mb — 20 Mq for our standard Salpeter IMF and 
this simple analysis would then predict the correct order 
of magnitude for the core collapse time (tcc/irh(0) ~ 0.1; 
see Figure 121). Note, however, that this simple dynami- 
cal friction picture corresponds to one or a few massive 
objects traveling through a uniform background of much 
lighter particles, a description that provides, at best, only 
a rough approximation of the real situation under con- 
sideration here (see Appendix). 

1.3. Goals of this Study 

In this paper, the first of a series, we consider a wide va- 
riety of initial cluster models and we investigate the evo- 
lution of all systems until the onset of core collapse. Re- 
sults from numerical simulations that include stellar col- 
lisions and track the growth of a massive object through 
successive mergers of massive stars during core collapse 
will be presented in a subsequent paper. This work is 
in progress and pre liminary results h ave already been 
reported elsewhere llRasio eTaOllOOl. Here we deter- 
mine how the core collapse time is related to various ini- 
tial parameters including the IMF (Sec. 3.2 and 4.2) and 
the central concentration and tidal radius of the cluster 
(Sec. 3.1 and 4.3), and we also include the possibility of 
initial mass segregation (Sec. 3.3 and 4.4). We then make 
a comparison with the stellar evolution time and derive 
limits on cluster initial conditions to allow the develop- 
ment of runaway collisions and the possible subsequent 
formation of a central massive BH (Sec. 5). We also de- 
rive from our calculations an estimate of the total mass 
of the stars that participate in the runaway collisions, 
which provides an upper limit to the final BH mass. 

We do not include stellar evolution in our calculations 
since our aim here is to investigate dynamical processes 
taking place before even the most massive main-sequence 
stars in the cluster have evolved. We do, however, study 
the effects of mass loss from stellar winds and their de- 
endence on metallicity (Sec. l4.5|l . In a subsequent paper 
Giirkan fc R,asioir2003() . we will study the evolution of 
"post-collapse" star clusters in which a central IMBH is 
assumed to have formed early on (through the runaway 
collision process). In particular, we will study the possi- 
bility that mass loss from the stellar evolution of the re- 
maining massive stars (those that have escaped the cen- 
tral collapse and runaway) could disrupt the cluster (cf. 
poshi ct al. 2001), thereby producing a "naked" IMBH. 
This is motivated by observations of ultraluminous X- 
ray sources in regions of active star formation (e.g., in 
merging galaxies) containing many young star clusters, 
but with the X-ray sources found predominantly outside 
of those clusters (Zezas & Fabbiano 2002). 

An important factor that could affect significantly the 
dynamics of core collapse in y oung star clusters is the 
presence of primordial binaries l|Fregeau et al.ll2003ar or 
the dynamical formation of binaries through three-body 
interactions (Gicrsz 1998; Hut 1985). As pointed out 
by Inagaki (1985), the formation rate of hard three- 
body binaries w ould be accel erated in clusters with a 
mass spectrum l|Heggiel[T975l) . In this first paper we 
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do not take into account the presence of binaries in 
the cluster. This is justified because we do not ex- 
pect binaries to play an important role until after the 
onset of core collapse, when the central density in- 
creases suddenly. However, collisions should become 
dominant immediately after the onset of core collapse. 
Our expectation, based on several previous theoretical 
studies of physical collisions during interactions of hard 
primordial binaries and for three-body binary forma- 
tion, is that the presence of binaries in the core will 
in fact increase c ollision rates , thereby helping to trig- 
ger the runaway tBacon ^t al.l lT996: Ch ernoff fc Huaiiel 
119961 iFregeau et al.ll2003bD . This is also supported by 
the results of direct A^-body simulations showing that, 
in smaller systems containing N > 10^ single stars, col- 
lisions indeed occur predominantly through the inter- 
actions of three-body binaries form ed at core collapse 
(jPortegies Zwart fc McMillanll200^ . 

2. NUMERICAL METHODS AND SUMMARY OF 
PREVIOUS WORK 

Numerical methods for investigating the dynamical 
evolution of star clusters include direct iV-body integra- 
tion, solutions of the Fokker-Planck equation by direct 
(finite-difference) or Monte Carlo methods, and gaseous 
models (for a review, see [Heggic & Hut 2003). Here 
we refer to direct iV-body integrations simply as 'W- 
body simulations," direct integrations of the Fokker- 
Planck equation as "Fokker-Planck (FP) simulations," 
and Fokker-Planck simulations based on Monte Carlo 
techniques as "Monte Carlo (MC) simulations." Note, 
however, that our MC simulations, in which the cluster 
is modeled on a star-by-star basis, are in fact another 
type of "A^-body simulations." Each approach offers dif- 
ferent advantages and disadvantages for understanding 
core collapse and massive black hole formation in star 
clusters with a realistic mass spectrum. 

2.1. Summary of Previous Numerical Work 

iPortegies Zwart fc McMillanI l|2002D carried out N- 
body simulations starting from a variety of initial con- 
ditions for clusters containing up to ~ 6 x 10^ stars. 
They found that runaway collisions driven by the most 
massive stars can happen in sufficiently dense clusters. 
Their results apply directly to small star clusters con- 
taining ~ 10"* — 10^ stars. However, in such a small 
cluster, any realistic IMF typically contains only a very 
small number of massive stars. For example, a Salpeter 
IMF with minimum mass TOmin = 0.2 Mq and maximum 
mass mniax = 120 M0 contains a fraction ~ 3 x 10~^ 
of its stars above 60 Mq. The dynamical role played by 
massive stars can therefore depend strongly on the to- 
tal number of stars in the cluster. In addition, in small 
systems, the dynamical evolution might be dominated 
by the random behavior or the initial conditions of just 
a few very massive stars. Consequently, to investigate 
runaway collisions in larger systems such as super star 
clusters or proto-galactic nuclei, realistically large num- 
bers of stars must be used in numerical simulations. The 
computational time required for direct integration of an 
A^-body system over one crossing time tdyn scales as 
(N^-'^^ on parallel machines, if o ne can adjust the num- 
ber o f CPUs to N optimally; see lSpurzem fc Baumeardj 
|2003). Since the relaxation time is ~ {N/ In N) ^dyn (see 



Sec. II. 2|) . the scahng of the total CPU time of A^-body 
simulations is nearly as steep as A'^'^. Currently, using a 
state-of-the-art CRAPE-6 board to accelerate the com- 
putations (Makino 2001, 200^), the evolution of a cluster 
containing 10^ stars with mmax/^min — 1000 and no pri- 
mordial binaries can be integrated up to core collapse in 
about one day (H. Baumgardt, private communication). 
However, including primordial binaries or increasing the 
number of stars would still lead to prohibitively long com- 
putation times, especially for a parameter study where a 
large number of integrations are required. 

A wide mass spectrum also leads to increased compu- 
tation times for FP simulations and gaseous models^. In 
these methods, a continuous mass spectrum is approx- 
imated by discrete mass bins. In the gaseous model, 
the most time-consuming operation of the algorithm is 
the inversion of large matrix whose dimension is pro- 
portional to the number of equations, itself proportional 
to A^comp- Hence the computing time increases like^ 
Tcpu oc A'^^omp- For FP codes, there is one diffusion 
term coupling each component to all others, leading to 
Tqpv oc A^comp- These steep scahngs limit A'^comp to at 
most 20 to avoid computation times longer than a few 
days. 

Even with a small number of mass bins, FP simu- 
lations have yielded important qualitative results. In- 
agaki and collaborat ors investigate d two- and multi- 
component systems (Tnagaki '1985^ llnagaki fc Saslawl 
^985; Inagaki & Wivanto 1984). They found that the en- 
ergy transfer within the core is an important process for 
determining the onse t of core collapse. FP simul ations 
by Inagaki (1983) and lChernoff fc WeinbergI 1(199(1) show 
that a mass function with mmax/^min — 10 — 15 needs 
to be discretized into about 15 components or more to 
obtain an accurate value of the core collapse time. A 
coarse r discretization leads to an a rtificially slow evolu- 
tion. iChernoff fc Weinberg) l)1990|) successfully tracked 
the energy transfer between different components and 
demon strated the importance of this process. Quinl^ 
(|1996l[) used FP simulations to follow the evolution of 
single- and two-component systems. He used a unique 
setting for his two-component clusters, with each com- 
ponent having a different initial spatial distribution. His 
aim was to study the interaction of galactic nucleus con- 
taining mainly dark (compact) objects with the bulge of 
the galaxy. In one model he considered a structure with 
"inverse initial mass segregation," where the central nu- 
cleus is made of objects much lighter than the normal 
stars composing the bulge. He showed that this situa- 
tion also leads to highly accelerated evolution as the more 
massive stars get trapped by the nucleus throu gh dynam- 
ical fr iction and undergo rapid core collapse. iTakahashU 
l(1997|) also investigated the evolution of clusters with a 
mass spectrum using FP simulations. His simulations 
were two-dimensional (in phase space) and therefore he 
was able to study the development of velocity anisotropy. 

We have been informed (H. Cohn & B. Murphy 2004, private 
communication) that latest isotropic FP codes require only a few 
seconds of computation to reach core collapse, for as many as 20 
mass components. 

^ Note that, in principle, one could split each step into separate 
"Poisson" and "Fokker-Planck" parts, in a way similar to what 
is done in FP codes, hence reducing the cost to Tcpu o<: Af^^^^p 
(R. Spurzem, private communication). 
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Gaseous models have the advantage of being very fast 
(for A'^comp < 20) but they include the greatest number 
of simplifying assumptions and so require independent 
checking. The good agreement shown with our MC sim- 
ulation for a Plummer model with Salpeter IMF in Fig- 
ure El is encouraging. Note, however, that the gaseous 
model calculation shown in Figure |21 used 50 mass com- 
ponents. This number of components is exceptionally 
high, leading to a total computation time in excess of 
three weeks, much longer than most MC runs. It was 
chosen to be able to follow core collapse to a very ad- 
vanced stage, until the most massive component dom- 
inates the collapsing core. We note that the evolution 
of the gaseous model to core collapse is faster by some 
30 % and exhibits a more gradual contraction of the in- 
ner region. The reasons for these small discrepancies will 
be investigated in future work, where additional com- 
parisons between Monte Carlo a nd gaseous model ca l- 
culations will also be presented ijFreitag et al.l I2003aj) . 
For the time being, suffice it to mention that for multi- 
mass clusters, there are basically two adjustable (di- 
mensionless) parameters in the gaseous model, one set- 
ting the effective thermal conductivity and the other the 
timescale for energy exchange between components (A 
and Acq , see 'Giersz & Suurzcm' "1994| iLouis fc SpurzernI 
[iS91; Spurzem & Takahashi 1995). For the model plot- 
ted in Figure we used the standard values of these 
parameters (A — 0.4977 and Acq = 1) which were es- 
tablished for clusters with a different structure and mass 
spectrum and some adjustment (preferentially through 
comparisons with A^-body runs) may be required. 

A direct comparison between FP and gaseous mo dels 
was also carried out bv lSpurzem fc Takahashil l)1995D . for 
two-component clusters, and also resulted in good agree- 
ment. Recently, a hybrid code has been developed by 
l^icrsz & Spurzem (2003), combining a gaseous model 
with MC techniques. In this approach the single stars 
are represented by the gaseous model, while primordial 
binaries are followed with a Monte Carlo treatment. A 
similar hybrid treatment could be applied to the prob- 
lem we are studying here, but with massive (single and 
binary) stars included in the Monte Carlo component, 
and lower-mass stars represented by a gaseous model. 

2.2. The Monte Carlo Code 

The solution of the Fokker-Planck equation with an 
orbit-averaged Monte Carlo method is an ideal compro- 
mise for the problem at hand. It can handle a suitably 
large number of stars and a wide mass spectrum can be 
implemented with very little additional difficulty. Most 
importantly, as in direct A^-body integrations, MC sim- 
ulations can implement a star-by-star description of the 
cluster. This allows the inclusion of many important pro- 
cesses such as collisions, binary interactions (including 
primordial binaries) , stellar evolution (and the accompa- 
nying mass loss) , as well as the effects of a massive central 
object with relative ease and much higher realism com- 
pared to direct FP simulations or gaseous models. We 
will incorporate the effects of all these additional pro- 
cesses on cluster evolution in the subsequent papers of 
this series. 

The MC code we have used to obt ain the main results 
of this paper is described i n detai l bvlJoshi et alJ l)2000D . 
It is based on the ideas of lHenonI l)1971alblll973j) and in 



m any respects it i s very simil a r to MC codes deve loped 
bvlStodoikiewicdl llT98a IT981 . IGierszl lll998ll2fl0l . and 
iFreitag fc'^en3 " l|2001f) . In the rest of this section, we 
give a brief summary of the numerical method. 

The main simplifying assumption is the Fokker-Planck 
approximation, in which relaxation processes are as- 
sumed to be dominated by small distant encounters 
rather tha n strong encounters with large deflections 
(jBinnev fc Tremainc 1987: Soitzcr 1987). The dynam- 
ical evolution can then be treated as a diffusion process 
in phase space. Following the individual interactions be- 
tween the stars, as in direct A'^-body simulations, is com- 
putationally expensive. However, the average cumula- 
tive effect on a star in a given amount of time can be 
characterized by diffusion coefficients in phase space. To 
compute the relaxation, the timestep for the numerical 
evolution has to be chosen smaller than the relaxation 
time. Since the relaxation time is normally shortest at 
the center, we choose our timestep to be a fraction of the 
central relaxation time. This ensures that the relaxation 
is followed accurately throughout the cluster. 

Another important simplification in the MC method 
is the assumption of spherical symmetry. The position 
of the particles is represented by a single radial coordi- 
nate, r, and the velocity is represented by radial, Vr^ and 
tangential, Vt^ components. The specific angular momen- 
tum, A, and the specific energy, of a star with index 
i, are given by: 

A, = vt,^ r,; , E, = U{r,) + ^(w^ ,, + ^2^), (7) 

where U (r) is the gravitational potential at a given point. 
The assumption of spherical symmetry implies that the 
potential at all points can be computed in a time pro- 
portional to the number of particles A^, rather than N"^ . 

In Henon's algorithm, the evolution is simulated by re- 
producing the effect of the cluster on each star by a single 
effective encounter in each timestep. At every iteration 
the two integrals of the motion Ai and Ei characterizing 
the orbit of each star are perturbed in a way that is con- 
sisten t with the value of the diffusion coefficients l)Henonl 
To conserve energy, this perturbation is realized 
by a single effective scattering between two neighboring 
stars. The square of the scattering angle, /3, is propor- 
tional to the timestep chosen. We choose our timestep 
such that sin^(/3/2) < 0.05 at the center of the cluster. 
Choosing too large a timestep will lead to a saturation 
effect and the relaxation will proceed artificially slowly. 
Using too small a timestep, on the other hand, not only 
increases the computation time, but also can lead to spu- 
rious relaxation (see the Appendix for a discussion of 
spurious relaxation effects). 

After the perturbation of Ai and Ei, the stars are 
placed at random positions between their apocenters and 
pericenters using a probability distribution that is pro- 
portional to the time spent at a given location on their 
new orbit. For a star with index i, the apocenter and 
pericenter distances are calculated by finding the roots 
of 

2E, - 2U{r) - ^ = 0. (8) 

This random placement is justified by the assumption 
of dynamical equilibrium, i.e., the evolution of the sys- 
tem does not take place on the crossing (or dynamical) 
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timescale, but rather on the relaxation timescale. The 
only important point about assigning a specific position 
to a particle on its orbit is that its contribution to den- 
sity, potential, interaction rates, etc., has to be estimated 
correctly. After all stars have been placed at their new 
positions, the potential is recalculated and the whole cy- 
cle of perturbation is repeated. 

This method can be modified so that the timestep is 
a fra c tion of the lo cal relaxation time jFre itag & Bcri3 
ppoT; Hcnon 197^. IStodolkiewicl iITq^ 1198611 and 
[Gicrsz (1998, 2001D divided the system into zones re- 
sulting in an approach intermediate between a fixed and 
smoothly varying timestep. Dividing the system into 
radial zones also allows the MC method to be paral- 
lelized efficiently for use on multi-processor machines 
ijJoshi et al.ll2000(l . Another possible modification uses 
a scaling of the units such that each particle in the sim- 
ulation can represent an entire spheric al shell of many 
ident i cal stars rather than a single star ijPreitag fc Bend 
120011 12002t IHenonllT97Ta|) . It should also be noted that, 
although the effects of strong encounters between stars 
on relaxation are assumed to be negligible compared to 
weak, more distant encounters, they can be incorporated 
by estimating their rate of occurrence in a way similar to 
physical collisions ("Freitag & Benz 2002^ or interactions 
with binary stars (Fregeau et al. 2003a). 

Our MC code has been used previously to study many 
fundament al d ynamical processes such as the Spitzer 
instabilitv ([p atters et al. 2000) and mass segregation 
l|Fregeau et'a]l2002 ) for simple two-component systems, 
as well as the evolution of systems with a continu- 
ous but fairly nar row mass spectrum of evolving stars 
l|Joshi et alJl200lD . A difficulty introduced by a broad 
continuous IMF (with a large mmax/TOmin ratio) is the 
necessity of adjusting the timestep to treat correctly en- 
counters between stars of very different masses. When 
pairs of stars are selected to undergo an effective hyper- 
bolic encounter as described above, one has to make sure 
that the deffection angle remains small for both stars. 
In situations where the mass ratio of the pair can be 
extreme, one has to decrease the timestep accordingly 
llStod61kiewiczll982 ^. In practice, for the simulations de- 
scribed here, we find that the timestep has to be reduced 
by a factor of up to ~ 500 compared to what would be 
appropriate for a cluster of equal-mass stars. We discuss 
further the applicability of orbit-averaged MC methods 
to systems with a continuous mass spectrum and large 
"T-max/^min ratio in the Appendix. 

3. INITIAL CONDITIONS AND UNITS 

The characteristics of core collapse and the subsequent 
runaway collisions depend on the initial conditions for the 
cluster. These initial conditions include the total num- 
ber of stars, the IMF, the initial spatial distribution of 
the stars (density profile and, possibly, initial mass seg- 
regation), and the position of the cluster in the galaxy, 
which determines the tidal boundary. As we shall see, 
the most important initial parameters are the slope of 
the IMF, the maximum stellar mass, and the initial de- 
gree of central concentration of the cluster density profile. 
We have used a wide variety of initial conditions, both 
to test the robustness of our findings and to establish 
the dependence of our results on these parameters. As 
a typical reference model we use an isolated Plummer 



sphere with a Salpeter IMF and stellar masses ranging 
from TOmin — 0.2 Mq to minax = 120 AIq. We then ex- 
plore variations on this model by changing the initial 
cluster structure, the IMF or the number of stars. 

3.1. Density Profile 

We have exami ned three families of mode ls: Plummer 
and King models ijBinnev fc Tremainell987l Sections 2.2 
and 4.4) J which have a core-halo stru cture, and 7-models 
l)Dehnen|ll993l: ITremaine et al.llT99^ . which have a cusp 
near the center. All these models have a characteristic 
radius given by 



ap 



?2/3 



l)i/Vh ~ 0.766 rh. 



9a2 



47rG/9o' 



(9) 
(10) 



^7 = (2'/^'"''' - l)^h, (11) 
for Plummer, King, and 7-models respectively. In these 
formulae, rh is the half-mass radius, po is the cen- 
tral density, and cr is a King model parameter. We 
show the density profiles corresponding to these vari- 
ous models in Figure 0] Here Wq is the dimensionless 
central p otential, related to t he concentration parame- 
ter llBinnev fc Trem ainc '1987'. Fig. 4-10). Other useful 
quantities characterizing the various density profiles are 
given in Table The initial model used for each of our 
MC simulations is listed in the second column of Table|21 
A general procedure for p roducing these models is given 
bv lFreitag fc Benz ("2002) ; a less general but simp ler pro- 
cedure for the Plummer model is given bv iAarseth et alJ 

3.2. Initial Mass Function 

Since we expect the conditions leading to core collapse 
to depend sensitively on the IMF, we have carried out 
simulations with a wide range of IMF parameters. How- 
ever, it is generally established that, at least for the high- 
mass end of the spectrum, a universal IMF very close to 
the simple Salpeter-like power-law is obeyed. This is in- 
dicat ed bo th by ob servations an d by theoretical calcula- 
tions ijClark c 2003: Krour)al2002L and references therein). 

In our simulations we assign the masses of individual 
stars using a sampling procedure. We first choose a ran- 
dom number, X, from a uniform distribution between 
and 1. For a simple power-law IMF with dN oc m~" dm 
between TOmin and m^ax we calculate the corresponding 
mass using 



m{X) 



l + X 



1 



1/(1-") 



(12) 

where the value a = 2.35 would correspond to a Salpeter 
IMF. 

In addition to simple power la ws, f or two mod- 
els we have used t he Miller-Scalo ()1979(1 and Kroupa 
ijKroupa et alJl99'^ IMFs, which are steeper at the high- 
mass end of the spectrum and shallower at the low-mass 
end. Both of these IMFs can be represented by broken 
power laws. For the Miller-Scalo IMF, 

fm-i'' 0.1<m<1.0 

dN 



dm 



,-3.3 



1.0 <TO< 10.0, (13) 
10.0 < m 
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Table 1. Properties of clusters 



Cluster 


PO 




a 


n 






Mc 






Plummer 


1.167 


0.532 


0.589 


oo 


0.769 


0.417 


0.192 


0.093 


0.0437 


King Wo = 1 


0.454 


0.534 


1.517 


2.568 


0.858 


0.670 


0.321 


0.110 


0.1134 


2 


0.530 


0.526 


1.003 


2.800 


0.849 


0.612 


0.281 


0.108 


0.0930 


3 


0.652 


0.518 


0.749 


3.134 


0.839 


0.543 


0.238 


0.106 


0.0722 


4 


0.860 


0.510 


0.576 


3.625 


0.827 


0.465 


0.195 


0.104 


0.0523 


5 


1.252 


0.504 


0.438 


4.362 


0.814 


0.382 


0.1546 


0.101 


0.0348 


6 


2.112 


0.503 


0.320 


5.471 


0.804 


0.293 


0.1171 


0.100 


0.0205 


7 


4.526 


0.511 


0.2146 


6.987 


0.812 


0.2032 


0.0830 


0.101 


0.00997 


8 


13.742 


0.530 


0.1253 


8.344 


0.872 


0.1211 


0.0531 


0.112 


0.00368 


9 


55.671 


0.558 


0.0649 


8.374 


0.980 


0.0633 


0.0307 


0.134 


0.00106 


7 = 1.0 


oo 





1/3 


oo 


0.805 








0.100 





1.5 


oo 





1/2 


oo 


0.851 








0.108 






Note. — Definitions of quantities listed in this table: pq is the central mass density; ctq is the central one-dimensional velocity dispersion; a is 
the characteristic radius (see Eas. lolllH : rt is the tidal radius; is the half-mass radius; rc — [Qctq /(47rGpo)] is the core radius; Mc is the mass 
enclosed by r^; trh is the half-mass relaxation time (cq.|^; and f^c is the central relaxation time fca. llBi . All quantities arc given in A^'-body units, 
except that i^h and trc are given in "Fokker- Planck units" (see Section l3.4l . 
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Fig. 4. — Density profiles for various initial models used in our 
simulations. The Plummer model is shown by a solid line, 7-models 
corresponding to 7 = 1, 1.5 and 2 are shown by dotted lines, and 
King models with Wo = 3, 6, 9 and 12 are shown by dashed lines. 
The radius is in units of the characteristic length scale parameter of 
each model and the density profiles are normalized to unity at that 
radius. King models with lower Wo values have steeper profiles 
outside r = ax, while 7-models with higher 7 values have higher 
p/p{a) near the center. 



and for the Kroupa IMF, 




0.08 <m< 0.5 
0.5 < TO < 1.0 

1.0 < TO 



(14) 



where all numerical values are in solar mass. 

However, for the sake of computational conve- 
nience, we prefer to implement a somewhat different 
parametrization of these distributions. To generate the 
mass spectra corresponding to these IMFs we use the 



functions ijEggleton et al.lll989l: iKroupa et cd]ll993(l 

0.19X 



m{X) 



(1 -X)0-75 -1- 0.032(1 -X)0-25 



(15) 



and 



m{X) ^ 0.08 



0.19X 



1.55 



omx 



0.6 



{l-Xf 



(16) 



for the Miller-Scalo and Kroupa mass functions, respec- 
tively. In both cases if m{X) < TOmin or m{X) > rrimax 
the result is discarded and a new random number is gen- 
erated. The values in equation IKill correspond to the 
choice of ai = 1.3 in Table 10 of Kroupa e t al. ( 1 993j) . 

In most of our simulations we have used N = 1.25 x 10^ 
stars. For all mass functions, this implies the pres- 
ence of many massive stars, allowing us to resolve fully 
the higher end of the mass spectrum. For example, 
for niniin = 0.2 M0 and TO^ax = 120 Mq, with N = 
1.25 X 10'', we have > 60 stars with to > 100 M© for 
a Salpeter IMF. The results from our simulations are 
therefore not much affected by random fluctuations in a 
small number of very massive stars. Note, however, that 
in much smaller systems, containing perhaps only ^ 10** 
stars (as in the Arches and Quintuplet clusters near our 
Galactic center), these small number effects and random 
fluctuations may indeed play a dominant role in deter- 
mining the dynamical fate of the few most massive stars 
in the system. 

3.3. Initial Mass Segregation 

Initial mass segregation in star clusters (i.e., the ten- 
dency for more massive stars to be formed preferentially 
near the cluster center) is expected to r esult from star 
formation feedback in dense gas clouds ijMurrav fc LinI 
11996(1 or from competitive gas accretion onto proto- 
stars and mer gers b etween them (Bonncll & Bate 200!^ 
iBonnell et aLll200H) . There is also some observational 
evidence for initi al mass segregation in both open and 
globular clusters llBonnell fc Daviesll 19981 Ide Griis et alJ 



l2?)03tlRaboud fc Mermilliodlll998D . 

We have considered the possibility of initial mass seg- 
regation in a few of our MC simulations. We adopt a 
simple prescription whereby we increase the average stel- 



lar mass within a certain radius 



For r < 



rather 
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than sampling from an IMF with fixed m^i^, we ran- 
domly choose between two values of mmin, one that is 
used for the outer part of the cluster and another one 
that is larger. For r > r^g, we follow a similar procedure, 
this time changing mmax- The average mass within r^g 
is larger by a factor Cms and outside rms is smaller by a 
factor C^g, with respect to a cluster without initial mass 
segregation. We choose Cms and Cms such that the over- 
all average stellar mass in the cluster does not change 
with these modifications. Changing the average stellar 
mass within any region of the cluster would of course in 
general leave the system out of dynamical equilibrium. 
To maintain virial equilibrium, the mass density profile 
must also be preserved. We achieve this by modifying 
the number density of stars appropriately. 

Our initial conditions for models with initial mass 
segregation are summarized in Table |31 Here q is the 
initial cluster mass fraction within r^g- This implies 
Cms = (1 - - qCnis)- Increasing q or Cms rep- 

resents more extended or more pronounced initial mass 
segregation. 

Our prescription for initial mass segregation allows the 
formation of massive stars in the outer parts of the clus- 
ter, as well as the formation of lighter stars in the inner 
parts, but the more massive stars are more likely to be 
found near the c enter. This makes our ap proach differ- 
ent from that of iBonnell fc DaviesI l)1998l) . who put all 
massive stars closer to center. 

3.4. Units 

For all our nume rical calculatio ns, we adopt the stan- 
dard iV-body units l)Henonll971a|) : we set the initial total 
cluster mass M = 1, the gravitational constant C = 1, 
and the initial total energy Eq = —1/4. In the tabulation 
of the results we also use the initial half- mass relaxation 
time trh(O), given by equation as the unit of time for 
comparison with other work in the literature. 

The conversion to physical units is done by evaluat- 
ing the initial half-mass relaxation time in years. For 
example, for the Plummer model, we can write 

trh{0) ^ 330 Myr 

^ / A^/lQg \ (ap_y^^ f M \-^'^ 
\ \xi{-icN)/\n{l{)^)) \\vc) 

(17) 

Similar expressions can also be obtained for King mod- 
els and 7-models by use of the quantities in Tabled For 
a single-component model (containing equal-mass stars) 
the value of 7c in th e Coulomb logarithm can be cal- 
culated theoretically ijFarouki k, Salr)eterlll994j) . but for 
a system with a wide mass spectrum it must be de- 
termined by comparing to d irect A'^-body integrations. 
iGiersz fc Heggid l|1996l Il997j) carried out such compar- 
isons and found 7c — 0.015. Our own comparison with 
a recent A^-body result for > 10^ stars with a wide mass 
spectrum led us to adopt the value 7c = 0.01 (H. Baum- 
gardt, private communication). 

For processes occurring in the central parts of the clus- 
ter, the central relaxation time is a more relevant quan- 
tity, 

0-3^ 

^"^'^^ = 4.887rG2 ln(7ciV)n(m)2' ^^^^ 



where ctsd, n and (m) are the 3D velocity dispersion, 
numbe r density and average stellar mass at the cluster 
center ()Spitzeilll987l eq. 3.37). The iV-body unit system 
only specifies unambiguously dynamical times, which are 
independent of N . In this system, relaxation times are 
proportional to 7V/ln(7cA^). It is therefore useful to de- 
fine also the so-called "Fokker-Planck time unit," which 
is the iV-body time unit [t^yn = GM^/'^{-AE(^)-^/^) 
multiplied by iV/ln(7ciV). 

The most important physical properties of all our ini- 
tial cluster models, including trh(O) and trc(O), are given 
in Table □(in TV-body and FP units). 

4. RESULTS 

The initial conditions and main results of all our MC 
simulations are summarized in Tables |21 and |3I All mod- 
els have initially N — 1.25 x 10^ stars, except for Mod- 
els 2r, s, b, c, and d, which have varying N between 
3 X 10^ and lO'^. The maximum value of A'' is set in prac- 
tice by the available computer memory and N 10*" 
corresponds to 2 GB of available memory for our code. 
The run for Model 2d took about two weeks of CPU time 
to complete on a 2.8 GHz Pentium 4 Linux workstation. 
More typical runs for iV = 1.25 x 10^ took around 20-30 
CPU hours. The close agreement between the outcomes 
of Models 2r - 2d confirms the expectation that our re- 
sults should be independent of the number of stars in the 
system, at least for sufficiently large N to avoid small- 
number effects in the cluster core. 

4.1. Mass Segregation and Core Collapse 

As expected, all models with a Salpeter-like IMF and a 
wide mass spectrum undergo core collapse considerably 
faster than any single-component cluster (cf. Sec. I1.2|l . 
In Figure |S1 we show the evolution of various Lagrange 
radii, as well as the average stellar mass inside these 
radii, for our reference model (Model 2; same as shown 
in Figure In contrast to the evolution of a single- 
component model, the inner Lagrange radii remain al- 
most constant until the very abrupt onset of core col- 
lapse at i/irh(0) — 0.07. A more detailed view of core 
collapse is shown in Figure|Sl Here the time axis has been 
replaced by the central potential depth, which increases 
monotonically in time and provides a natural stretch near 
core collapse. Note that the core collapse time for a par- 
ticular run can be determined very accurately, to within 
^ 0.01%. However, random fluctuations in the realiza- 
tion of each initial condition for a particular model lead 
to a much larger physical "error bar" on ice (see below). 

The rapid increase of the average stellar mass inside the 
innermost Lagrange radii seen in Figures and El is an 
indication of significant mass segregation. In fact, very 
significant mass segregation takes place throughout the 
evolution of the system. This can be seen in the upper 
panel of FigureEl which shows the evolution of half-mass 
radii for stars in v arious mass bins (compare with Fig. 1 
of S pitzer fc Sh uU 1975). It is clear that the rate of mass 
segregation in each mass bin, measured by the slope ams 
of the corresponding half-mass radius r{t) in that bin, is 
very nearly constant from t — Q all the way to core col- 
lapse (the least-square straight-line fits, constrained to 
r/rh(0) = 1 at i = 0, arc shown in the figure). This rate 
can be positive or negative. For Model 2 we find that 
stars more massive than about 5 Mq drift inward on av- 
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Table 2. Initial conditions and main results 





Initial Structure 


IMF 






(Mq) 


(m) 


tcc/trh(,U; 


iWcc/iWtot 


55 


Plummer 




PL, 


a 


— 


■ 


0.2-120 


6.57 


0.287 




50 


Plummer 




PL, 


a 




— 1.7 


0.2-120 


2.74 


0.131 


0.0024 


1 


Plummer 




PL, 


a 




—2.0 


0.2-120 


1.28 


0.0899 


0.0018 


51 


Plummer 




PL, 


a 




—2.2 


0.2-120 


0.87 


0.0700 


0.0020 


2r 


Plummer 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.0716'' 




2s 


Plummer 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.0702 




2 


Plummer 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.0700*" 


0.0020'" 


2b 


Pluuimer 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.0700*" 


0.0019*" 


2c 


P lummer 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.0706 


0.0020 


2d 


l-* In TTl W1 f^V 
JT ILllillliCl 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.0720 


0.0018 


52 


Plummer 




PL, 


a 




-2.5 


0.2-120 


0.58 


0.0719 


0.0022 


3 


1-^ In TT1 m tfiv 

JT ILllillliCl 




PL, 


a 




-2.7 


0.2-120 


0.48 


0.0696 


0.0018 


53 


P lummer 




PL, 


a 




-3.0 


0.2-120 


0.40 


0.0834 


0.0012 


4 


Plummer 




Kroupa 




0.2-120 


0.96 


0.0858 


0.0014 


5 


T-*! 1 1 TT1 m 

1 ILllliiliCl 




Miller-Scalo 


0.2-120 


0.71 


0.0723 


0.0022 


28 


P lummer 




PL, 


a 




—2.35 


0.2-360 


0.72 


0.0795 


0.0020 


27 


Plumm.er 




PL, 


a 




-2.35 


0.2-240 


0.71 


0.0760 


0.0020 


24 


In TT1 m 

IT ILllliiliCl 




PL, 


a 




-2.35 


0.2-90 


0.68 


0.0664 


0.0014 


20 


P 111 TTl in 

IT 1 LliillliCl 




PL, 


a 




-2.35 


0.2-60 


0.67 


0.0786 


0.0024 


21 


P 111 TTl in 
JT ILliiiiliCl 




PL, 


a 




-2.35 


0.2-20 


0.62 


0.156 


0.0028 


22 


Plumm.er 




PL, 


a 




-2.35 


0.2-8 


0.56 


0.478 


0.0010 


25 


t--' 111 TTl in 

IT 1 LliillliCl 




PL, 


a 




-2.35 


0.2-5 


0.53 


0.805 


0.0016 


23 


111 TTl m 
IT 1 LliillliCl 




PL, 


a 




-2.35 


0.2-2 


0.45 


2.20 




26 


t-' 111 TT1 m CiV 
iT ILliliiiiCl 




PL, 


a 




-2.35 


0.2-1 


0.37 


4.29 




30 


xviii^, ' ' — 1 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.152 


0.0020 


36 


-LVillg, VV — -i 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.151 


0.0014 


35 


TCino- Wrx — 9 
ri.iiig, I'V'o — -^1 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.134 


0.0020 


37 


ilvilig, KkQ — ^ 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.129 


0.0014 


11 


j:\.iiig, VV Q — o, 


fit 


PL, 


a 




-2.35 


0.2-120 


0.69 


0.107 


0.0022 


10 


King, Wo = 3 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.110 


0.0016 


31 


King, Wo = 4, 


(i) 


PL, 


a 




-2.35 


0.2-120 


0.69 


0.0779 


0.0022 


38 


King, Wo = 4 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.0778 


0.0022 


32 


King, Wo = 5, 


(i) 


PL, 


a 




-2.35 


0.2-120 


0.69 


0.0561 


0.0024 


oy 


King, Wo = 5 




PL, 


a 




-2.35 


u.z— izu 


u.Dy 


U.Uozu 


U.UUzu 


12 


King, Wo = 6, 


(i) 


PL, 


a 




-2.35 


0.2-120 


0.69 


0.0336 


0.0014 


40 


King, Wo = 6 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.0322 


0.0020 


33 


King, Wo = 7, 


(i) 


PL, 


a 




-2.35 


0.2-120 


0.69 


0.0163 


0.0014 


41 


King, Wo = 7 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.0150 


0.0018 


34 


King, Wo = 8, 


(i) 


PL, 


a 




-2.35 


0.2-120 


0.69 


0.00545 


0.0010 


42 


King, Wo = 8 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.00577 


0.0012 


13 


King, Wo = 9, 


(i) 


PL, 


a 




-2.35 


0.2-120 


0.69 


0.00135 


0.0010 


43 


King, Wo = 9 




PL, 


a 




-2.35 


0.2-120 


0.69 


0.00138 


0.0012 


6 


7=1, Hernquist 


PL, 


a 




-2.35 


0.2-120 


0.69 


< 10-4 




7 


7 = 1.5 




PL, 


a 




-2.35 


0.2-120 


0.69 


< 2 X 10"'' 





Note. — All models have N = 1.25 x lO'' stars except for Model 2r (N = 3 X 10^), Model 2s (Af = 6 X 10^), Model 2b (N = 2.5 X lO''), Model 2c 
{N — ^ X 10^), and Model 2d {N — 10^). Isolated King models are indicated by (i). When an entry is missing in the last column (•■■), we were 
not able to determine Mcc reliably for that model. 

''These results are obtained by averaging over 20 runs. 

*"These results are obtained by averaging over 10 runs. 



erage, while less massive stars drift outward. Even for 
the most massive stars, we do not find that the mass seg- 
regation rate is proportional to the average mass in the 
bin (in contrast to what would be expected for a tracer 
population of massive st ars driven by s imple dynamical 
friction; cf. Sec. 1.2 and lFreeeau et al] l2002'). Instead, 
the following simple expression provides a good fit (to 
within a few percent) to the observed mass-dependence 
for Model 2: 

amsirh(O) = ao exp(-m/TOf) -|- 02, (19) 

where m is the average stellar mass and the best-fit pa- 
rameters are ap = 9.45, nif = 21.9 Mq, and a2 = —8.07. 
Note that the mass segregation rate actually approaches 
a constant for large m (but of course it is unphysical to 
extrapolate beyond rrimax)- In the lower panel of Fig- 



ure \7\ we show the average mass within rh , rh /2 and 
rh/4. The steady increase of the average mass in each 
region is further indication that the mass segregation not 
only starts immediately, but also continues until core col- 
lapse. For a smaller number of stars, the average mass 
wi thin a given radius can reach saturation (cf. Figure 2 
of lBonnell fc Davie£lll99^ . Our results for larger N do 
not show this saturation. 

Initially and throughout the evolution until core col- 
lapse, the cluster maintains a core-halo structure. How- 
ever, when the heaviest stars start dominating the core 
and the Spitzer instability occurs, this initial structure 
is lost. We demonstrate this behavior in Figure|Sl where 
the evolution of the gravitational potential profile is 
shown. We have checked that the final structure and, in 
particular, the formation of a cusp are not dominated by 
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Table 3. Initial conditions and results of simulations with initial mass 
segregation. 



0.04 0.06 

t/t,,(o) 

Fig. 5. — Evolution of the Lagrange radii (bottom) and the av- 
erage mass within Lagrange radii (top) for Model 2. The Lagrange 
radii are given in units of the initial half-mass radius and time in 
units of the initial half-mass relaxation time. A more detailed view 
concentrating on core collapse is given in Figure 151 
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Fig. 6. — Same as Figure 1^ but concentrating on the evolution 
of the cluster near core collapse. The horizontal axis now gives 
the value of the central potential, normalized to its initial value. 
The vertical dashed line marks the core collapse time, <cc/trh(0) = 
0.068. Dotted lines on either side indicate a ±0.01% change in 
this quantity, illustrating how precisely the onset of core collapse 
can be determined for a particular run. The corresponding mass 
fraction in the collapsing core is A/cc/A/tot = 0.0018 for this run. 



Model 


q 


fms/rrh(0) 


C'ms 


; 

^ms 


tcc/trh(0) 


Mcc/Mtot 


mOl 


0.3 


0.69 


1 


2 


1.094 


n 0588 


0.0018 


ml7 


0.3 


0.69 


1 


5 


1.273 




n 0090 


m02 


0.3 


0.69 


1 


8 


1.522 


0.0443 


n 0099 


ml8 


0.3 


0.69 


2 


1 


1.892 




0022 


m04 


2 


0.55 


1 


2 


1.053 


0.0637 


n 00*^0 

u. uuou 


ml5 


0.2 


0.55 


1 


5 


1.143 


0.0512 


0.0026 


m05 


0.2 


0.55 


1 


8 


1.250 


0.0498 


0.0020 


ml6 


0.2 


0.55 


2 


1 


1.379 


0.0439 


0.0022 


m06 


0.2 


0.55 


2 


4 


1.538 


0.0399 


0.0022 


m08 


0.1 


0.40 


1 


2 


1.023 


0.0664 


0.0016 


ml3 


0.1 


0.40 


1 


5 


1.059 


0.0588 


0.0018 


m09 


0.1 


0.40 


1 


8 


1.098 


0.0558 


0.0016 


ml4 


0.1 


0.40 


2 


1 


1.139 


0.0560 


0.0020 


mlO 


0.1 


0.40 


2 


4 


1.184 


0.0519 


0.0018 


mil 


0.1 


0.40 


3 





1.286 


0.0506 


0.0022 


ml2 


0.1 


0.40 


3 


6 


1.406 


0.0471 


0.0016 



Note. — All initial models arc isolated Plummer spheres containing 
1.25 X 10^ stars. Here q is the mass fraction initially contained within 
^ms- Inside this radius, the average stellar mass is larger by a factor 
Cms compared to Model 2. Outside this radius the average stellar mass 
is smaller by a factor C^^. See Sec. 3.3 for details. 




0.04 

t/t,,(o) 

Fig. 7. — Mass segregation in Model 2. In the upper panel we 
show the evolution of the half-mass radii of stars in various mass 
bins. From top to bottom, the limits of the mass bins, in solar 
mass, are 0.2, 0.32, 0.55, 0.94, 1.62, 2.77, 4.7, 8.1, 14, 30, 50, and 
120. In the lower panel we show the average stellar mass within, 
from top to bottom, rij/4, ^1^/2, and rj,. 



small number effects but are instead the result of many 
massive stars participating in the core collapse. We illus- 
trate this in Figure IHl where close to 200 innermost stars, 
which constitute 0.2% of the total mass, are shown ex- 
plicitly. 

Another important point, which can be deduced from 
the evolution of the Lagrange radii, as well as from Fig- 
ure |H1 (see the lower two curves), is that the outer parts 
of the cluster are not much affected by mass segregation 
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Fig. 8. — Evolution of the gravitational potential for Model 2. 
The potential U is given in units of the absolute value of the initial 
central potential. The initial profile is shown by the dashed line. 
The profiles just before core collapse and right after core collapse 
are shown by dotted and solid lines, respectively. The positions of 
the stars in the collapsing core are marked individually. For this 
particular run, the core collapse took place at t = 0.069 trh(O). 
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0.1 - 



core MF 

Salpeter IMF 



10 
M/Mg 



100 



Fig. 9. — Comparison between the mass spectrum of stars in 
the collapsing core (solid line) and the Salpeter IMF (dashed line) 
for Model 2. Numbers for the core mass function are obtained by 
averaging over three successive snapshots. The Salpeter IMF is 
normalized so that both curves account for the same total number 
of stars. 



and core collapse, except for the disappearance of the 
most massive stars. 

As expected, the stars in the collapsing core have a 
mass distribution much richer in heavy stars than the 
IMF. We compare the mass function at core collapse 
with the IMF, for m > IMq, in Figure El The mass 
function shown in this figure is obtained by choosing the 
innermost 200 stars at three different times all before and 
within 5 x 10~^trh(0) of core collapse time. 

Mass segregation continues until the core is dominated 
by massive stars. At this point, the massive stars near 
the center dynamically decouple from the rest of the clus- 
ter and go into collapse as a separate subsystem. The 
stars outside this subsystem act as an energy sink, and, 
as a result of their energy gain, expand away from the 
center of the cluster. In a single-component system, as 
the collapse proceeds, the subset of stars that partici- 
pate in the collapse becomes smaller and smaller. Here 
instead the stars that participate in the collapse at the 
edge of the core are heavier than the surrounding stars 
so they can more efficiently give their energy away and 
remain together. As a result, in systems with a wide 
mass spectrum and a steep enough IMF, there is a clear 
separation between stars in the collapsing core and stars 
outside the collapsing core, almost suggesting a real con- 
densation process. The total mass in the collapsing stars 
is a crucial property of these systems, clearly represent- 
ing an upper limit to the mass of any BH that could 
form eventually through the gravitational collapse of a 
runaway merger remnant. 

To estimate the total mass in the collapsing core, Mcc, 
we proceed as follows. At every timestep we calculate 



and record the Lagrange radii for various mass fractions. 
For each Lagrange radius, this provides typically a few 
thousand data points per run. The innermost radii ex- 
hibit a great amount of noise for the lowest mass frac- 
tions. We remove this noise in two steps. First, we take 
arithmetic averages over 60 points and reduce the total 
number of points to 50 — 100. Then for each point and 
using the closest 6 points we fit a cubic polynomial us- 
ing least squares and we evaluate this polynomial at the 
corresponding time. This way we both smooth the La- 
grange radii further and estimate their derivatives. When 
the system is going into core collapse, a decrease in the 
collapse rate is an indication that the results provided 
by our code are beginning to be dominated by numerical 
errors. We stop the simulation when this happens and 
define the core collapse time as the last point before this 
behavior is observed. We then find the innermost La- 
grange radius that has a positive derivative at the time 
of core collapse and estimate the mass of the stars that 
participate in core collapse as the mass enclosed by this 
radius. An investigation of Figure by eye verifies that 
this method produces reasonable results. In practice, we 
track many more Lagrange radii than shown in this fig- 
ure, at intervals of 0.02% between 0.08% and 0.36%. 

Because of small-number fluctuations (and, in par- 
ticular, the intrinsic noise in the innermost Lagrange 
radii) there is a statistical uncertainty in all the numbers 
quoted in Table|21 This uncertainty is not unphysical, as 
real systems will be affected similarly by fluctuations in 
the number and specific properties of a relatively small 
number of massive stars. To test the robustness of our re- 
sults and estimate this uncertainty we have repeated 10- 
20 simulations with different random seeds for Model 2 
and some of its variants with different numbers of stars 
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(Models 2r, 2s, and 2b). For these models the numbers 
in Table El are obtained by averaging the results of the 
various (physically equivalent) runs. The standard devi- 
ations obtained for the values of fcc/irh(0) and Mcc/Mtot 
for our typical Model 2 are about 5% and 20%, respec- 
tively. The larger uncertainty in Mcc/Mtot is a result 
of the high sensitivity of this quantity to noise in the 
innermost Lagrange radii. 

We cannot study the evolution of the cluster past core 
collapse without treating in detail the dynamics of the 
central stars during collapse, which is beyond the scope 
of this paper. In future work we will include a de- 
tailed treatrnent of s tellar collisions in the core (Paperll, 
iFreitag et alJl2003bl) . Alternatively, one could use a hy- 
brid method and treat the central part of the cluster 
with a direct A^-body approach (Light man fc McMillanI 
Il985j) . However, the simplest approximation may be to 
introduce an effective boundary condition at some very 
small radius, ^todolkicwicz ( 1983), in his model B2, took 
this approach, but his implementation was not conserv- 
ing energy. 

4.2. Dependence on the IMF 

The sharp onset of core collapse for our reference 
model, shown in Figure |S1 is a result of the Spitzer insta- 
bility. This instability is driven by the segregation of the 
heaviest stars toward the center and their dynamical de- 
coupling from the rest of the system. As indicated by the 
results of Section 14.11 the ratio of maximum to average 
stellar mass in the IMF, mmax/ (™): is an important pa- 
rameter setting the timescale for the onset of instability. 
There are various ways to study the dependence of our 
results on this parameter. One can use an IMF different 
from Salpeter, e.g., Miller-Scalo or Kroupa. In addition, 
when using a power-law IMF, changing the slope a or 
the maximum mass mmax will obviously alter the value 
of the ratio Wmax/ {m) . 

Our results from simulations for a large number of 
models exploring these various alternative IMF are pre- 
sented in Figure ^] Note that for mmax/ (m) < 40 
our results suggest a relation tec oc (mms.^/ ( rn))~^'^, 
also o beye d by all Fokker-P lanck models from llnaeakil 
()1985|) and lTakahashil l)1997l) with a < 3.5 and a num- 
ber of mass components sufficient to ensure proper sam- 
pling of the IMF. However our computations extend to 
much higher values of TOmax/(w) than those works. Be- 
yond rrimax/ (w) ~ 50, a domain reached by any real- 
istic IMF, the core collapse time approaches a constant 
~ 0.15 ire (0). Therefore our main conclusions appear to 
be independent of the details of the IMF as long as the 
number of massive stars in the system is large enough. 

For small values of rrimax/ i'm), not only the timescale 
but also the very nature of the collapse changes. In these 
systems the evolution timescale, i.e., the relaxation time, 
of the subsystem that can decouple from the rest of the 
cluster is no longer small enough that the evolution of 
lighter stars can be neglected. The Lagrange radii for 
these models behave similarly to the single-component 
case, i.e., there is no clear separation between collaps- 
ing and expanding Lagrange radii. Consequently, we do 
not give values of AfccM^tot for these models in Table |21 
iQuinla n (1996), using FP simulations, also found that, 
for a clear decoupling, the relaxation time of the sub- 
system has to be short compared to that of the other 
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Fig. 10. — Dependence of the core collapse time on the shape 
and width of the IMF. Initial Plummer models were assumed in all 
cases. The horizontal axis shows the ratio mmax/ (m) of maximum 
to mean stellar mass in the IMF. The square is for our standard 
model with a Salpeter IMF, ?72min = 0.2 Mq , and TTimax — 120 M0. 
The cross and the plus are for the Kroupa and Miller-Scalo IMF, 
respectively, with the same limits. The triangles are for power- 
law IMF with various exponents (-1.4, -1.7, -2.0, -2.2, -2.5, -2.7, 
-3.0) and the same limits. The circles are for Salpeter-like IMF 
with varying upper limits (mmax = 1, 2, 5, 8, 20, 60, 90, 240, 
and 360 Mq). For comparison, FP results from|lnagaki 11985., his 
table 2) are plotted with small solid dots. Models with the same 
IMF slope (a) and A^comp = 5 (number of discrete mass com- 
ponents in the FP simulations) are linked together with dashed 
lines. For each value of a, a model with a relatively large mass 



contrast (mmax/m^ 



10) was computed with A^c< 



15. 



Its core-collapse time is shorter than a corresponding model with 



5 because the mass function is better represented. These 



models are linked to the other points with dotted lines. The 
three multi-mass models (A^comp = 10) from iTakahashi. 1.199^ are 
also plotted with small open dots (joined with dash-dotted line). 
They correspond to a = 3.5, 2.5, and 1.5 (in order of increasing 
mmax/ (m)). 



stars. 

4.3. Dependence on Initial Cluster Concentration 

One expects the core collapse time to depend strongly 
on the initial density profile of the cluster, and, in par- 
ticular, on the central concentration. A simple and sys- 
tematic way to examine this dependence is to use differ- 
ent Ki ng mode ls with varying concentration parameter 
or Wo ijBinnev fc Tremainei ll987. Fig. 4-10). For single- 
component clust ers this has been done using a variety of 
method s ( Eins ellT99l I Joshi et alJl200ll iKim et alJl2002t 
'Quinlan< ll99(i^) . 

We have carried out a number of simulations using 
King models with Wq = 1 — 9 as initial structure. Our 
simulations include both isolated clusters (with no tidal 
boundary enforced, even though the initial models are 
truncated) and clusters with a tidal boundary (assuming 
a circular orbit in a spherical galactic potential). This 
tidal boundary is initially chosen to be at the tidal radius 
of the King model and then adjus t ed as the cluster loses 
mass llChernofr fc Weinberg|[T99?i I.Toshi et alJl2nnU) . In 
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Fig. 11. — Comparison of core collapse times for various King models with a broad IMF and with single-mass stars, (a) The ratio of core 
collaps e time to ha l f-m ass relaxation time. For the single-component mo dels (top curves), we include results obtained with Fokker-Planck 
codes fEinsel"199a; "Kmi et al. 20021: lOuinlan.i,1996 'l. Monte-Carlo codes IJoshi et ani20011 Freitag, unpublished) and one Wo = 3 A^-body 
simulation by Baumgardt (2001). Solid and open symbols indicate models with and without tida l truncation, respectively, (b) The ratio 
of core collapse time to central relaxation time. For the sake of clarity, only data from ljoshi et all I2001f) are plotted for single-component 
clusters. 



Figure [TTk we plot the ratio of core collapse time to ini- 
tial half-mass relaxation time for our simulations, along 
with results for single-component clusters. The ratio 
icc/irh(0) is much smaller for systems with mass spec- 
trum, as would be expected. We also see that the core 
collapse time for these systems is independent of the pres- 
ence of a tidal boundary, as expected when core collapse 
is driven by local processes within the core rather than by 
global energy transfer. This idea is strongly supported by 
the results in Figure lTTb . which shows the comparison of 
the ratio of core collapse time to central relaxation time, 
for our simulations and single-component models. With 
respect to single-component models, the ratio icc/irc(0) 
for systems with a broad mass spectrum shows very little 
variation. 

We have also run two simulations for 7-models, us- 
ing 7 = 1.0 (Hernquist model), and 7 — 1.5. Initially, 
models with 7 < 2 have vanishing central velocity dis- 
persion and, for 7 < 2, they exhibit a central "tem- 
perature inversion" ijTremaine et al.l|l994|) . For single- 
component clusters, the central region at first undergoes 
rapid gravothermal expansion until it b ecomes isother - 
mal and "normal" core collapse can start l|Quinlanll99fiD . 
Models with 7 < 2 also have zero initial central relax- 
ation time^. Our results for King models would then 
suggest that core collapse induced by central mass seg- 
regation should be extremely fast, if it were not for the 
opposite effect of the temperature inversion. For 7 = 2, 
mass segregation seems to have the upper hand in the 
competition and it proves impossible to resolve core col- 
lapse because extreme mass segregation appears nearly 

® Models with 2.0 < 7 < 3.0 have infinite frc( O) because their 
veloc ity dispersion rises like 1/r near the center JTremaine et al] 



instantaneously during the MC runs. Models with 7 = 1 
and 7 = 1.5 expand at first and then evolve to core col- 
lapse very rapidly. The core collapse time and collapsing 
core mass are very hard to determine numerically for 
these models. In Table|21 we give only an upper limit on 
their core collapse times. Obviously, these values are so 
short that their physical relevance is unclear. It is hard 
to imagine through which process such a cluster could 
be created if one wants it to be virialized but without 
initial mass segregation because these conditions impose 
constraints on the formation timescale unlikely to be sat- 
isfied in real systems. In addition, note that the local 
dynamical time in a 7 = 1 model approaches a constant 
non-zero value near the center, while, formally, the re- 
laxation time goes to zero there. This suggests that such 
models can only provide an approximate description of 
real clusters, where finite-number effects will play a key 
role near the center, so that the question of determin- 
ing the evolution in the limit of very large N is not well 
posed. 

4.4. Initial Mass Segregation 

Naturally we expect that any initial mass segregation 
in a cluster should lead to an even shorter core collapse 
time as the heavier stars are starting their life closer to 
the center on average, and therefore do not need as much 
time to concentrate there through mass segregation. 

Our results from MC simulations with initial mass seg- 
regation (Sec. 3.3) are presented in TableOland FigureEl 
As expected, we find that tcc/irh(0) decreases with in- 
creasing C,„s (i.e., when stars become more and more 
massive on average in the inner region), and with in- 
creasing q (i.e., when the increase in average stellar mass 
takes place over a larger central region). However, the 
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Fig. 12. — Dependence of the core collapse time tcc/irh(0) on 
initial mass segregation parameters Cms and q (see Table IS! . 



values of Mcc/Mtot do not appear to be affected signif- 
icantly (within expected random fluctuations from run 
to run). The implication for BH formation through run- 
away collisions is that the final BH mass may not be 
affected by initial mass segregation, while the condition 
for runaway growth to occur (see Sec. 5) will be relaxed 
as the core collapse takes place even earlier. 

4.5. Effects of Stellar Evolution 

The simulations described so far consider only the dy- 
namical evolution of the cluster, neglecting the stellar 
evolution entirely. We have also carried out a number of 
simulations to understand how stellar evolution and the 
accompanying mass loss can modify the dynamical evo- 
lution. The st ellar evolution trea tment we have adopted 
is that of Belc zvnski et al.l i20n2l) . based on the approx- 
imations of H urlev et alJ l)200Cl() . 

Stellar evolution introduces a new physical clock in the 
system, independent of relaxation. Therefore it is nec- 
essary to specify the physical scale of an initial cluster 
model (e.g., the half-mass radius in pc) before starting 
a simulation so that the relaxation time irh(O) can be 
calculated in years and related to the stellar evolution 
timescale. The only other parameter to be specified is 
the metallicity Z, which plays an important role in the 
calculation of stellar evolution mass loss. In our treat- 
ment of stellar evolution, all wind mass loss rates are 
proportional to Z'^^'^. Since stellar evolution introduces 
two new parameters in our initial models, a full explo- 
ration of the initial parameter space is clearly impossible. 
Fortunately, we will see that the effects of mass loss on 
the evolution to core collapse are rather unimportant, 
and so a systematic study is unnecessary at this point. 

We first illustrate how post-main-sequence evolution of 
massive stars can prevent core collapse. In Figure [T^ we 
show the evolution of a system similar to Model 2 with 



tj.h(0) = 60Myr and Z = 10^^. At such low metaflic- 
ity, there is little mass loss on the main sequence, so one 
expects tec - 0.07 x 60 Myr = 4.2 Myr without evolu- 
tion beyond the main sequence. However, massive stars 
evolve off the main sequence after about 3 Myr, before 
core collapse has occurred. Once they evolve off the 
main sequence, they typically lose up to 2/3 of their 
mass rapidly. At the end of their life, which extends 
about 10 % beyond their main-sequence lifetime, these 
stars undergo supernova explosions and their cores col- 
lapse to BH. In the upper panel of Figure [T51 we plot the 
number of these BH in the cluster. The significant mass 
loss during the late stages of stellar evolution causes the 
cluster core to expand so that collapse is prevented and 
the system then goes into long-term dynamical evolution. 

Varying Z from 10""* to 0.02 (the solar value), we 
found that this reversal of core contraction occurs for 
any metallicity when massive stars are allowed to evolve 
off the main sequence. Higher metallicities of course lead 
to increased mass loss and an even stronger tendency for 
the cluster core to expand. We have not considered the 
case of Z = 0, i.e., a cluster of Pop III stars. These 
stars may collapse to BH t hat incorporate ess entially the 
entire initial stellar mass iHeger et al J 120031) . The dy- 
namical evolution of the cluster would then be unaffected 
by mass loss. Whenever a runaway is avoided, a dense 
cluster of relatively massive primordial BH would then 
form near the center. This is an intriguing possibility 
to keep in mind for future consideration, although the 
IMF of Pop III stars is essentially unknown (but see, 
e.g.,^bia ct al. 2001; Nakamura & Umcmura 2001) and 
recent hydrodynamic simulations suggest that the first 
stars may actuall y form isolated rather than in clusters 
llAbel et alJl2n0l . 

We have also examined the possibility that mass loss 
on the main sequence may already be a source of indirect 
heating strong enough to reverse core collapse or, at least, 
delay it significantly. Our results indicate that wind mass 
loss from main-sequence stars alone can never prevent 
core collapse and that the delay introduced remains very 
small, even for high metallicities. It is of course always 
possible to fine-tune the relaxation time so that the mass 
loss slows down the evolution just enough to allow a few 
massive stars to evolve off the main sequence and stop 
contraction. But such models can only represent a very 
small domain in the parameter space of initial conditions. 
For example, in a system like Model 2 with irh(O) = 
40 Myr, the core collapse time increases only by about 
5% for Z = 0.02. For t,.h(0) = 50 Myr, the evolution 
is just delayed long enough to prevent core collapse for 
Z — 0.02, but not for Z = 0.001. One has to increase the 
relaxation time to irh(O) = 60 Myr for stellar evolution 
to prevent core collapse for all metallicities. 

In conclusion, we find that post-main-sequence mass 
loss can be strong enough to prevent early core collapse, 
but that it does not significantly tighten the condition 
(on irh(O)) for core collapse to occur. Mass loss from 
main-sequence stars always plays a negligible role. 

5. SUMMARY AND DISCUSSION 

It has long been realized that a mass spectrum ac- 
celerates the dynamical evolution of dens e star clusters 
throug^h the process of mass segregation l)AarsetUll966l: 

is a consequence of 



througn the p rocess ot mass segreg; 
IHenoniirOTlat IWielenlllQTl . This 
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Fig. 13. — Prevention of core collapse due to stellar evolution 
mass loss. The lower panel shows the evolution of the Lagrange 
radii. The upper panel gives the number of black holes in the clus- 
ter (formed by evolving massive stars). The time is given in units 
of trh(O) for the bottom and middle axes, and in Myr along the 
top axis. The first massive stars leave the main sequence around 
3 Myr. 



the statistical tendency of two-body gravitational en- 
counters to establish energy equipartition between stars 
of different masses. For any realistic IMF, equiparti- 
tion can never be achieved because the massive stars 
quickly form a separate self-gravitating subcluster (a sys- 
tem with negative effective heat capacity) as they seg- 
regate near the center by tra nsferring kineti c energy to 
lighter stars ilnagaki fc Sasb- w 1985; Inagaki fc Wivantol 
1 1 984HSpitzeil 1 9691 1 Vishniad 1978 1. This instability leads 
to a rapid core collapse: after a finite time ice, the cen- 
tral density of stars would actually become infinite in 
the absence of finite size effects (both physical collisions 
between stars and the finite number of stars in the core). 

With the exceptio n of a fe w rec ent TV- 
body sim ulatio ns l lPortegies Zwa rt et al.l Il999t 
iPortegies Zwart fcTicMillanl |2003|) . all previous in- 
vestigations of core collapse in clusters with a mass 
spectrum have considered a relatively narrow range of 
stellar masses (TOmax/'TT-min ^ 100). This restriction 
is appropriate for old globular clusters in which stellar 
evolution probably had time to remove all massive stars 
before significant relaxation took place. In contrast, 
in the present study, we ignored stellar evolution to 
concentrate on systems in which, by assumption, core 
collapse occurs before even the most massive stars (with 
TO ~ 100 M©) leave the MS. 

Through a set of high-resolution Monte Carlo sim- 
ulations of clusters with a variety of IMF and struc- 
tural parameters (concentration, presence or absence of 
tidal truncation or initial mass segregation), we estab- 
lished two important results. The first concerns the time 
needed to go deep into core collapse {tec)- As the evo- 
lution to core collapse is driven by relaxation processes, 



^cc can always be written as proportional to the initial 
half-mass or central relaxation time (trh(O), trc(O)). The 
proportionality constant does not depend on the size of 
the cluster or the number of stars, as long as they are 
numerous enough to avoid small- number effects. For iso- 
lated clusters with an IMF of realistic slope (a ~ 2.5 in 
the high-mass range), we find 

icc/U(0) oc ^i-^ with n = and 5 ~ 1.3, (20) 

(to) 

as /i increases from 1 (single- mass) to 50. This 
result extends and agr ees nicely with previous wor k 
based on FP simulations l|Inagakiil 985t iTakahashi 1 997f) . 
From simple arguments about dynamical friction, one 
would naively expect a l inear relation with 5 = 1 
l)Binnev fc Tremaind Il987f ). The steeper dependence 
found in our numerical results may be related to the fact 
that the instability causing the core collapse is triggered 
only after some critical number, N„, of sufficiently mas- 
sive stars have drifted to the core by dynamical friction. 
If one has to go out to radius i?cr to find this number 
of heavy stars, the required timescale will be of order 
l^~^t-cci{Rci)- For systems with increasing ^, it is reason- 
able to expect iVcr and, consequently, i?cr and iroi(-Rcr) 
to decrease. While this provides a plausible explanation, 
further investigations will be needed to confirm it, in par- 
ticular stud ies making use of th e gaseous model of cluster 
dynamics iFreitag et al.ll2003aj) . 
For /i > 50, the core collapse time saturates to 

ice 0.15 t,c(0), (21) 

a key result given that any realistic IMF is likely to have 
/X > 100. This value of tec is at least two orders of 
magnitude shorter than the core collapse time for single- 
component clusters. Moreover, we find this tcc/irc(0) 
ratio to hold for all the models we have considered that 
have a finite irc(O). This includes King models (with 
or without tidal truncation) with Wq ranging from 1 to 
9, a sequence along which tcc/irh(0) decreases by more 
than 100. Our result that the core collapse time is 
set fundamentally by trc(O) rather than trh(O) refiects 
an important difference between single-component clus- 
ters and systems with a broad mass spectrum: in sys- 
tems with a mass spectrum, core collapse is driven ulti- 
mately by energy transfer occurring locally in the core; 
in contrast, relaxation in a single-component system is 
a global process taking place on all scales, resulting in a 
timescale comparable to irh (Inagaki 1985). For Plum- 
mer models or King models with moderate concentra- 
tion (Wo < 4), we find tec ^ 0.07 - 0.15trh(0), in good 
agreement with the results of A'^- body simulations by 
IPortegies Zwart fc McMillanI l)2002() '^. 

The main-sequence lifetime of very massive stars ap- 
proaches a constant value, of about 3 — 4 Myr, with only 
weak dependence on metallicity, rota tion, or the mass to 
of the star, as long as m > 50Mr7^ lIMaed er fcMe vnetl 
l200lt iMe^et fc Maededl200fi ISchaller et alJll993n t 
exhibits little variation with to because such massive 
stars are nearly Eddington-limited, hence L oc m, where 
L is the luminosity, and oc /cto/L ~ constant if the 

^ Although the fundamental relation is between tec and trc(O), 
trh is, by far, easier to estimate observationally. 
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fractional mass of the convective core, /c, does not de- 
pend too strongly on the stellar mass®. Therefore, the 
necessary condition for core collapse not being stopped 
by stellar evolution can be written 

trh(O) < 30Myr. (22) 

Interestingly, we also find that core collapse cannot be 
prevented by mass loss through stellar winds on the MS. 
Such mass loss can only delay the evolution to core col- 
lapse very slightly. Even for high mctallicities and hence 
strong winds, only clusters for which the core collapse 
time would otherwise be very close to the critical value 
(i.e., just below t^,) can see their fate changed through 
wind mass loss if ice increases to become slightly longer 
than t^. These clusters must represent a very small do- 
main in the parameter space of initial conditions. 

Our second important finding is that, as the collapse 
proceeds, the mass contained in the ever-shrinking core 
converges to a non-zero value, 

Mcc ~ 0.001 - 0.002 Mtot, (23) 

which depends only weakly on the properties of the initial 
cluster. This contrasts again with single-component clus- 
ters which, exhibiting self-similar collapse, have Mcc — > 

OjS i ~" ^ ^cc ' 

This value of Med Mtot remarkable agreement 

with the fractional mass of a dark object, presumably 
a massive BH, discovered at the center of galaxies (in 
which case, Aftot stands for the mass in the spheroidal 
component) and, possibly, globular clusters like M15 and 
Gl. This suggests that the collapsing core may well be 
the progenitor of a central IMBH. However, for a galac- 
tic nucleus, the requirement that tec — 0.15irc(0) < 
3 — 4Myr would imply that, either it evolved from a 
cluster with short (central) relax ation time to i ts pre sent- 
day state with > 10® yr ijLauer et al.l Il998j) . or 
that the growth of the central BH was predominantly 
from accretion of smaller star clusters with short re- 
laxation times (each one carrying an IMBH at its cen - 
ter: lEbisu^kiet al.y 200lV Hansen fc Milosavlievid2003j) . 
M15 and Gl also have long half-m ass relaxation times 
(~ 1 Gyr for M15.IDull et al.llT997l and 30 - SO Gyr for 
Gl. lBaumgardt et al.l2003bHMevlan et al.l2001(l indicat- 
ing that the process we envision could only have taken 
place if these clusters were born with much more cen- 
trally concentrated density profiles than observed today^ 
(for M15, an initial value of Wq > 8 seems required; see 
Fig. Ill() . However, there is no doubt that young star 
clusters with half-mass relaxation times shorter than the 
critical value of 30Myr do exist. The Arches cluster near 
the center of our Galaxy is the best establi shed exam- 
ple, with probably shorter than lOMyr ijFiger et al.l 

* If the mass segregation time scale is oc 1/m is compared to 
the "usual" relation i, oc 1/m'', one could conclude erroneously 
that massiv e stars never play an important role in core collapse 
iApplegatelllfiiSfill) . However, the approximate Xjiv? scaling of the 



Il999aj) . while the Quintup let cluster appears to He very 
close to the critical value l|Figer et al.lll999b(l . It is also 
possible that some super star clusters, the birthplaces of 
most stars in starburst en vironments, have sufficiently 
short relaxation times fe.g. IHo fc Filippen ko 1996). 

Although the final fate of the stars that participate in 
core collapse is not completely certain, very high rates 
of physical collisions are expected in the core. Primor- 
dial binaries are probably unable to stop the collapse and 
prevent these collisions. In smaller systems like globular 
clusters interactions with hard binaries may in fact in- 
crease the collision rate. In more extreme environments 
like (proto-)galactic nuclei, the velocity dispersion is so 
high (~ 10^ — lO^kms^^) that most primordial bina- 
ries are soft and will have been disrupted before they 
can play a role in the collapse. The main uncertainty 
is then whether these collisions, occurring at relatively 
high velocity, allow a massive star to gain mass in a run- 
away fashion, as suggested by analytical models (based 
on the dependence of the cross section on the mass and 
negle cting coUisional mass-loss; .Malyshkin fc Goodmaij 
120011 and references therein), or, on the contrary, grind 
it down progressively. Our preliminary results from MC 
simulat ions including a r ealistic treatment of stellar col- 
lisions ijRasio et al.ll2003l) indicate that the formation of 
a very massive star through runaway collisions is a likely 
outcome. More extensive calculations, as well as a de- 
tailed discussion of the likely fate of the merger remnant, 
will be presented in Paper II. 
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stellar lifetime applies only for m < 10 Mq. 

^ It is not clear that the present-day value of t^ji for a cluster, 
which can be estimated from observations, must be close to ti.h(0). 



APPENDIX 

APPLICABILITY OF THE MONTE CARLO METHOD 

The Monte Carlo method has been applied previously to two-component systems (see, e.g., iWatters et al.l 120001) 
and to systems with a continuous mass function covering a relatively narrow range (with m,„ax/'7imin ^ 30; see the 
references in Section However, one may question the applicability of the method to systems with a much broader 
mass spectrum. In particular, the nature of the energy transfer in these systems is different than in systems with 
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Fig. A14. — Test of spurious relaxation in Model 2. 
Here we show the evolution of the 10%, 50%, and 90% 
Lagrange radii, the central potential, the average mass 
within Lagrange radii, and the central density. All quan- 
tities are divided by the average of their first ten values 
in the run. 



small mmax/wniin ratio (See the discussion in Sec. 1.2). In addition, the gravitational potential is assumed to change 
smoothly within the cluster in the MC method. However, in a system containing a broad IMF, a relatively small 
number of very massive stars could introduce significant discrete changes in the slope of the potential. As the positions 
of these steps change randomly from one iteration to the next, one may worry that this could introduce a significant 
amount of spurious relaxation in the system. In this Appendix, we provide numerical tests and theoretical arguments 
showing that the MC method is indeed applicable to systems with a very broad mass spectrum, such as those under 
consideration in this paper. 



Spurious Relaxation 

In MC simulations, the orbits of the stars are calculated using a potential assumed to be spherical and smooth. 
This potential is constructed numerically from a radial distribution of point masses and therefore includes random 
fluctuations that can lead to spurious relaxation effects (Hcnon 1971a). To test the importance of these effects, we 
carried out simulations in which we turned off physical relaxation (by skipping the perturbation of stars). Ideally, in 
this case, the cluster should maintain a constant state indefinitely and none of its properties should change over time. 
In practice, however, it is sufficient to make sure that the changes are not significant compared to changes occurring 
in the presence of physical relaxation. 

In one of our tests, we let a MC simulation, for our typical Model 2, run for 3600 iteration steps, which is twice the 
number of steps normally required for this model to reach core collapse. We plot various relevant quantities showing 
the evolution of the system in Figure IA14I It is clear that the deviations in these quantities are much smaller than 
the changes due to physical relaxation, so we conclude that the effects of spurious relaxation is not significant. 

Relaxation and Dynamical Friction 

The energy transfer from massive to lighter stars and the segregation of massive stars to the center of the cluster 
are processes similar to sim ple dynamical friction for a heavy test particle embedded in a background of lighter stars 
(^inncv & Trcmain'3 ll987l Sec. 7.1). Some of the processes that govern the dynamics of dense stellar systems cannot 
be treated simply as part of two-body relaxation. Well-known examples include collisions, binary interactions, and 
three-body formation of binaries. A natural question to ask is whether dynamical friction is a process that can be 
modeled correctly by our implementation of two-body relaxation in the MC method. The answer is central to the 
applicability of the MC method to our problem. 

The frictional drag on a massive object is proportional to the mass d ensity of light background stars, and is inde- 
pendent of the masses of individual stars fsee lBinnev fc TremaindFigSTl eq. 7-18). So the effect of dynamical friction 
would be unaltered, for example, if every light star were replaced by two stars of half the mass. However, the treatment 
of relaxation in the MC method (Sec. 2.1) does depend on the mass ratio of the two particles considered and would 
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Fig. A15. — Orbital decay of a massive particle 
through dynamical friction. We plot the distance from 
the center, in A'^-body units, as a function of time in 
Fokker-Planck units (Sec. 3.4). In these units, the evo- 
lution does not depend on the value of 7c. We use our 
MC code to follow the decay of one massive object in 
a cluster of much lighter stars. The initial con ditions 
are identical to those used by ISpinnato et alJ 12002.) . 
The solid line shows the MC simulation result while the 
dashed line is the analytic solution from the standard 
Chandrasekhar theory of dynamical friction. See text 
for details. 



be altered by such a replacement. The ratio of relaxation to crossing time is given by 

tdyn In N 

which would evidently change if the stars in a system were replaced by many more stars with smaller masses. This 
might suggest t hat the MC m ethod cannot follow dynamical friction correctly. However, note that the relaxation time 
can be written l)Spitzedll987l eq. 3-37) 

3 

tr OC (A2) 
and the timescale for dynamical friction is proportional to ijBinnev fc TremaindllQS'^ eq. 7-18) 

tuF cx (A3) 
nm M 

Here, v is the average velocity of background stars, of mass m and number density n, and vm is the velocity of the 
massive test particle, of mass M . When vm is comparable to v (which is true, for example, if the massive particle is 
on a circular orbit; see below), we conclude that 

in V 

tnv oc — tr oc . (A4) 

The first proportionality implies that dynamical friction can be treated as a relaxation process. In the second, we 
show explicitly the dependence on the total mass density p of the background stars. This result i s also established by 
several numerical calculations (Bonnell & Davies 'lOOS"; Fre eeau et a l. 2002: Spit zer fc Shul]|ll97'5|l . 

To demonstrate that the MC algorithm is indeed able to handle simple dynamical friction, we have simulated 
a system consisting of a single massive object (mass 7712) in a cluster of much lighter particles (mass Wi , n umber 
iVi). The initial conditions are identical to those used for the Ni = 400 000 model of 'Spin nato et all l)2002|) . with 
m2/mi — 211, hence the total number of particles used for this test is 400 001. Spinnato e t alJ were interested in 
the spiral-in of an IMBH in a galactic nucleus. The massive particle is initially on a circular orbit at a relatively 
large distance from the center, r(0). For the power-law density profile used to represent the cl uster, p(r) oc r~'^ , one 
can predict anal yticall y the decay R{t) using Chandrasekhar's formula for dynamical friction (jSpinnato et al.ir2002L 
eq. 5). In Figure lA 151 we compare the spiral- in of the massive particle with the analytical solution. In the formula, 
we have set 7 = 2.3 to fit (visually) the density profile of the cluster. The initial iV-body set-up was provided by 
ISpinnato et al.l and converted directly into MC particles. The agreement between the MC run and the analytical 
formula is very satisfactory, especially considering that the density profile does not follow the exact power law and 
in fact evolves slightly during the simulation as a result of relaxation between light particles. Note that this good 
agreement does not validate the assumptions used in deriving Chandrasekhar's formula (e.g., neglecting large-angle 
scatterings, correlations between light particles, and gradients in the properties of the cluster) because the MC method 
relies basically on the same set of assumptions. In particular, 7c is a free parameter in both approaches. What this 
result demonstrates is that these assumptions are correctly implemented in our MC code. 
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